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Summary 


In  this  research,  the  s-version  finite  element  method  (s-FEM)  is  used  to  carry  out 
analyses  on  the  particulate  composite  materials  undergoing  progressive  damage.  Matrix 
material  is  assumed  to  be  a  polymeric  material.  Hard  particles  are  embedded  in  matrix 
material.  S-FEM  simplifies  the  modeling  procedures  for  particulate  composite  materials 
because  it  allows  us  to  build  finite  element  models  for  the  structure  as  whole  and  for 
particles  and  their  vicinities,  separately.  They  are  called  “global”  and  “local”  finite 
element  models.  When  particulate  composite  materials  are  modeled,  the  local  finite 
element  model  contains  a  particle  and  its  immediate  surrounding  region.  The  local  finite 
element  models  are  superposed  on  the  global  model.  By  adopting  the  s-FEM,  placing 
particles  in  matrix  material  became  a  trivial  task. 

Matrix  is  considered  to  suffer  from  a  material  damage  due  to  the  growth  and  nucleation 
of  microvoids.  Also,  the  composite  experiences  damage  due  to  particle-matrix 
dewetting.  The  former  is  accounted  for  by  the  use  of  a  continuum  damage  constitutive 
law.  The  later  is  by  a  cohesive  zone  model.  Two  kinds  of  continuum  damage  models  are 
used  in  this  research.  They  are  an  isotropic  and  a  separate  dilatational/deviatoric  damage 
constitutive  law  that  accounts  for  the  influences  of  hydrostatic  and  deviatoric  stresses 
separately.  The  constitutive  and  the  cohesive  zone  models  were  implemented  in  the 
s-FEM  computer  program. 

Numerical  analyses  were  carried  out  to  reveal  the  characteristic  deformation  behavior  of 
particulate  composite  such  as  the  influences  of  hard  particles  to  matrix  damage,  the 
influences  of  particle-matrix  dewetting  to  matrix  damage,  etc. 

The  results  revealed  that  the  separate  dilatational/deviatoric  damage  constitutive  model 
with  a  small  contribution  from  the  deviatoric  stresses  was  the  most  appropriate  model 
among  the  models  which  were  tested  in  this  study.  When  the  cohesive  zone  model  is 
assumed  at  the  interface  between  the  particles  and  matrix  material,  the  strength  of  bond 
of  the  cohesive  zone  determines  the  damage  mode  that  dominates  the  other.  When  the 
bond  is  strong,  the  matrix  damage  is  the  major  damage  mode.  When  the  bond  is  weak, 
the  dewetting  is  the  major  damage  mode. 

The  outcomes  of  present  research  revealed  some  characteristic  behavior  of  progressive 
damage  in  polymeric  particulate  composites. 
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1.  S-version  Finite  Element  Method  (S-FEM)  for  the  analysis  of 
particulate  composite  materials 


1.1  Equation  Formulations 

In  this  section,  we  discuss  about  general  formulations  of  the  s-version  finite  element 
method  (s-FEM).  It  is  assumed  that  particles  or  fibers  are  distributed  in  the  domain  of 
analysis  as  shown  in  Figure  1 .  Though  Figure  1  implies  that  the  particles  or  voids  are 
spherical  in  their  shapes,  there  is  no  such  restrictions  in  the  mathematical  formulation. 
The  second  phase  material  can  be  fibrous  or  any  others  in  their  shapes.  Though  the 
previous  formulations  of  the  s-FEM  (see  Fish  [1]  for  example)  assumes  only  one 
overlaid  model  to  be  superposed  on  the  global  model,  we  allow  any  numbers  of  finite 
element  models  to  be  superposed.  Then,  the  overlaid  models  are  allowed  to  overlap 
each  other,  as  depicted  in  Figure  2.  When  the  shapes  of  the  embedded  second  phase 
materials  are  the  same  or  similar  to  each  other,  the  same  local  finite  element  model  can 
be  used  repeatedly.  Therefore,  generating  a  model  for  the  composite  would  be  a  simple 
task. 

In  the  following  discussions,  the  regions  of  the  global  and  the  p- th  (p=  1,2,3,  1 1 1 L  M) 
local  finite  element  models  are  designated  to  be  C1G  and  CiLp ,  as  depicted  in  Figure  3. 
We  assume  that  there  are  a  total  of  M  local  model  regions.  The  displacements  are 
defined  based  on  the  shape  functions  of  elements  in  the  global  and  local  models, 
independently  [see  references  such  as  Bathe  [2]  and  Hughes  [3]  for  the  shape  functions 

of  finite  elements].  We  write  them  to  be  uf  =  uf(x)  in  and  ufp  =  ufp(x)  in  D.Lp , 

where  x  denotes  the  position  of  a  material  point.  At  a  point  which  is  not  inside  of  any 
local  model  regions,  the  displacements  m,  are  the  same  as  the  displacement  functions 
uf  of  aG . 


«/  =  u?  (*) 


(1) 


At  a  point  where  some  local  finite  element  models  overlap,  the  displacement  functions 
m,  are  given  by  the  sum  of  displacement  functions  of  the  overlapped  models.  For 
example,  at  a  point  where  the  local  models  nLp  and  Q.Lq  (1  <p,q<M,p*q)  overlap 
on  the  global  model  Q6 ,  the  displacements  ut  are  represented  by  the  sum  of  their 
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displacement  functions,  as: 


Uj  =  Wp  (x)+  ufp  (x)+  u^1  (x) 


(2) 


To  assure  the  continuities  of  displacements,  those  based  on  a  local  finite  element  model 
are  set  to  be  zero  at  its  outer  boundary.  We  let: 


u\p  =  0  at  dO.Lp 


(3) 


where  dCiLp  designates  the  outer  boundary  of  local  model  region  ClLp  . 

Stresses  at  a  point  are  written  in  terms  of  strains  through  Hooke’s  law  [see  Sokolnikoff 
[4],  for  example]. 


aij  ~  E ijkt£U 


(4) 


where  the  elastic  constants  Eijk,,  may  vary  within  the  solid  and  are  the  functions  of 
location  of  a  material  point. 


Eijkl  ~ 


EiM{x) 


(5) 


The  statement  of  principle  of  virtual  work  is  written  to  be: 

y^Em  l^dQ G  =  ]nG  +  Jan?  ^4d )  (6) 

j  ^ 

where  Su,  are  the  variations  of  displacements,  b,  are  the  body  force  per  unit  volume, 
tj  are  the  prescribed  traction  vector  on  the  traction  prescribed  boundary  80.^ .  The 
variations  of  displacements  Sut  are  assumed  in  the  same  manner  as  the  displacements 
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iii  >  by  the  superposition  at  material  points  where  they  overlap.  Thus,  Su,  are  written  to 
be  Sui  =  Suf (x)  where  no  local  models  overlap  on  the  global  model  and 

Sui=Suf{x)+Sufp{x)+Sufq(x)  where  local  models  ClLp  and  ClLq  (1<  p,q<M ,  p  *  q) 

overlap  on  the  global  model.  Suf  (l  <  p  <  m)  are  set  to  be  zero  at  the  boundary  dO.Lp 
of  nLp . 

Thus,  the  displacements  and  their  variations  are  substituted  in  the  statement  of  virtual 
work  principle.  After  some  algebraic  manipulations,  we  arrive  at: 


dSf 

dx  : 


du 


" ijki 


dx„ 


k  dQ° 


m  dSu 
p=l  ox 


du 


Lp 


“ ijki  ' 


dxf 


m'p  =  JqC  Sufbi<mG  +  \dnc  &,r,d(oQ,G )  (7) 


r  8Su‘P  F 


ijki 


8u<k  6D.Lp  +  Jj 


dXn 


dSufp 
dX  ■ 


~EW 


8ufP  J  M 
~^dO.Lp  +  I 

dXn 


9=1 

q*p 


Jq lp~^  " 


8Su\p 
dX  ■ 


~EijU 


r)uLq 

— < k—6D.Lp~Lq  =\0lp  SufpbjdQ.Lp 
dx ,  n 


(8) 

(p  =  1,2,3,---,m) 

From  the  left  hand  sides  of  equations  (7)  and  (8),  various  stiffness  matrices  are  obtained. 
Thus,  an  equation  can  be  written  in  a  matrix  form,  as: 


kg 

KC-L1 

G—L2 

kg~l  3 

kg~lm  " 

uG 

fg  ' 

ku-g 

Ku 

gLl-L2 

g  LI-L3 

g  Ll-LM 

uLX 

fLl 

KL  2-G 

g  L2-L1 

kl1 

g L2-L3 

£ L2-LM 

uL1 

fL2 

L3-G 

g  L3-LI 

g L3-L2 

kl 3 

K  L3~LM 

uL 3 

>  < 

pL3 

LM-G 

gLM-Ll 

j^LM-L2 

j^LM-L3 

...  Klm 

uLM 

flm 

Kg  is  the  ordinary  stiffness  matrix  for  the  global  finite  element  model  which  arises 
from  the  first  term  of  equation  (7)  and  KLp  (p  =  1,2,3,---,m)  are  those  for  the  local 
finite  element  models  arising  from  the  second  term  of  equation  (8).  kg~Lp 
(p  =  1,2,3, and  KLp~G  (p  =  1,2,3, are  the  coupling  stiffness  matrices  between 
the  global  and  local  finite  element  models,  arising  from  the  second  term  of  equation  (7) 
and  the  first  term  of  equation  (8),  respectively.  KLp~Lq  ( p,q  =  1,2,3 p*q)  are  the 
coupling  stiffness  matrices  between  the  local  finite  element  models  arising  from  the 
third  term  of  equation  (8).  FG  and  FLp  (p  =  1,2,3, -,m)  are  the  nodal  force  vectors 
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of  external  and  body  forces.  It  is  noted  that  the  matrices  have  the  properties  of 
klp-g  _  [Kc-Lp  J  anc|  KLP-Lq  _  ^ Lq-LpJ  Therefore,  the  coefficient  matrix  in  the  left 
hand  side  of  equation  (9)  is  symmetric. 


Unknown  nodal  displacements  are  obtained  by  solving  the  linear  simultaneous 
equations  (9),  and  the  displacement  field  in  the  composite  is  determined  (see  Okada,  Liu, 
Ninomiya,  Fukui  and  Kumazawa  [5]  for  the  full  details  of  the  solution  procedures). 


Figure  1  Schematic  views  of  composite  materials 


Figure  2  An  s-FEM  model  for  composite  material  in  which  each  fiber/particle 
and  its  immediate  vicinity  are  modeled  by  a  local  finite  element  mesh 
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Figure  3  Local  model  regions  (QLp(p  =  1,2,3,---,m))  which  are  superposed  on  the 
global  model  region  nG 

1.2  Numerical  Implementations  of  s-FEM  (Evaluation  of  stiffness  matrices) 

A  unique  feature  of  the  s-FEM  is  that  elements  in  global  and  local  finite  element  models 
overlap  each  other.  There  are  many  ways  for  elements  to  overlap,  as  shown  in  Figure  4. 
Elements  may  overlap  each  other  in  an  arbitrary  manner.  This  raises  serious  problems  in 
the  evaluation  of  stiffness  matrices  such  that  i)  when  a  coupling  stiffness  matrix  such  as 
K  Lp-Lq  js  formecj  elements  in  different  finite  element  models  may  partially  overlap 
each  other  (two-dimensional  illustration  is  presented  in  Figure  5)  and  ii)  more  than  one 
material  models  or  material  parameters  that  are  specified  by  the  overlapping  elements 
may  exist  at  a  point.  Some  special  care  must  be  given  to  overcome  these  issues.  In  this 
section,  how  we  can  overcome  these  issues  in  s-FEM  computer  implementation  is 
described. 
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(d) 


(e) 


(f) 


Figure  4  Many  ways  that  two  elements  overlap  each  other;  (a)~(c): 
Two-dimensional  examples  and  (d)~(f):  Three-dimensional  examples 


Figure  5  Overlapping  region  Q LPl  Lqi  of  the  regions  QLp‘  and  Q Lqi  of  two 

elements  (Lpi  and  Lpj) 

1.2.1  Evaluation  of  element  and  coupling  stiffness  matrix 

We  consider  a  typical  scenario.  Figure  5  illustrates  the  finite  elements  i  and  j  of  the  local 
mesh  regions  p  and  q.  They  are  designated  to  be  elements  Lpi  and  Lqj,  respectively. 
They  overlap  each  other  and,  therefore,  their  coupling  stiffness  matrix  is  formed.  The 
coupling  stiffness  matrix  [&  lp'~uu  J  can  be  written  to  be: 
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(10) 


k  Lpi-Lqj  ]  _  y  [BLpiJ 


~Lqj 


where  [b  Lpi  J  and  [b  Lqj  J  are  the  strain-displacement  matrix  for  the  elements,  and  [£] 
is  the  matrix  representing  the  elastic  constants  Eijk(  in  equation  (5).  QLpi~Lq>  is  the 

volume  of  overlapping  region.  When  the  coupling  stiffness  matrix  is  formed,  a 
numerical  integration  is  performed  for  the  overlapping  volume  nLp'~Lqj . 

However,  the  overlapped  region  nLpi~L®  may  have  a  complex  geometry.  It  is  almost 
impossible  to  explicitly  define  the  geometry  of  the  overlapping  region  and  to  apply  an 
ordinary  numerical  integral  scheme  such  as  Gasuss  quadrature.  In  this  study,  we 
perform  the  integral  based  on  one  of  the  overlapping  elements.  It  can  be  shown  to  be: 


k 


Lpi-Lqj 


a 


Lpi-Lqj  I 


(x)BLpi  T[E]BLqi\ mLpi 


(11) 


aLpi  Lqj(x)  is  a  scalar  function  whose  value  is  “1”  in  CiLpi  Lqj  and  “0”  outside  of 
Vt Lpi-Lqj  ,  as  indicated  in  Figure  5.  Numerical  integral  is  performed  based  on  element  Lpi. 
Therefore,  the  integrand  of  equation  (11)  has  a  sever  discontinuity  since  the  value  of 
aLP'-Lqi(x)  changes  abruptly  from  “1”  to  “0”  or  “0”  to  “1”.  Therefore,  an  ordinary 
Gauss  quadrature  is  unable  to  evaluate  the  integral  accurately.  In  order  to  circumvent 
this  problem,  we  developed  an  element  subdivision  scheme. 

The  element  subdivision  technique  is  illustrated  in  Figure  6  for  two-  and 
three-dimensional  problems.  As  shown  in  Figure  6  (a-i),  two-dimensional  elements  Lpi 
and  Lqj  partially  overlap  each  other,  and  the  integration  is  performed  based  on  Lpi.  First, 
element  Lpi  is  divided  into  4  sub-cells.  Then  each  subdivided-cell  is  checked  if  it 
intersects  with  the  edges  of  element  Lqj.  If  a  subdivided-cell  intersects  with  the  edges  of 
element  Lqj,  it  is  divided  into  4  sub-cells  again.  This  process  is  repeated  until  the 
smallest  sub-cell  becomes  small  enough  (typically  within  1%  of  the  volume  of  element 
Lpi).  The  processes  of  creating  sub-cells  are  shown  in  Figures  6  (a-ii)-(a-iv).  Then,  an 
ordinary  Gauss  quadrature  rule  is  applied  in  each  subdivided-cell.  The  same  approach  is 
adopted  in  three-dimensional  problems,  as  shown  in  Figure  6  (b-i)-(b-iv).  However,  it 
should  be  noted  that  proposed  methodology  is  not  computationally  efficient  and  takes 
much  longer  computational  time  than  the  ordinary  Gauss  quadrature,  because  many 
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integration  points  are  used  to  carry  out  the  numerical  integral. 


(b-i)  (b-ii)  (b-iii)  (b-iv) 

Figure  6  Element  subdivision  techniques  for  (a)  two-  and  (b)  three-  dimensional 
problems 


1.2.2  Material  constants  and  strain  histories 

In  an  ordinary  finite  element  method,  the  material  constants  are  assigned  to  finite 
elements  and  the  strain  history  information,  such  as  damage  paramters  are  stored  at  the 
integration  points.  In  present  s-FEM  analysis,  the  global  and  local  finite  element  models 
overlap  each  other  and  material  data  including  initial  strains  and  strain  history 
parameters  are  assigned  to  global  and  local  finite  element  models  independently.  This 
means  that  two  or  more  sets  of  material  parameters  exist  within  an  overlapping  region. 
First,  priority  orders  are  assigned  to  material  models  in  the  input  data.  Deformation 
history  data  is  stored  at  ordinary  integration  points  in  each  element  in  global  and  local 
finite  element  model.  When  two  or  more  finite  elements  overlap  at  a  point,  a  material 
model  having  the  highest  propriety  order  is  chosen  first.  If  there  were  two  or  more 
overlapping  elements  having  the  same  material  model  at  a  point,  deformation  history 
data  of  the  smallest  element  is  used  to  form  the  stiffness  matrices. 

For  example,  in  Figure  7,  a  two-dimensional  schematic  illustration  is  given.  There  are 
two  material  models  A  and  B  that  are  assigned  to  the  elements.  At  a  point  within  the 
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overlapping  region,  material  model  A  or  B  is  chosen  according  to  their  assigned  priority 
orders.  Even  when  an  ordinary  element  stiffness  matrix  is  formed,  abrupt  changes  in 
material  constant  may  occur.  The  element  subdivision  scheme  of  section  1.2.1  is  used 
for  such  cases. 

Next,  we  discuss  about  treatments  for  the  strain  history  parameters.  When  the  element 
subdivision  technique  is  adopted,  many  integration  points  just  for  numerical  integration 
are  generated.  If  one  tries  to  store  the  strain  history  data  at  each  one  of  them,  a  large 
amount  of  computer  memory  will  be  required.  Thus,  in  present  s-FEM  program,  each 
ordinary  integration  point  carries  the  strain  history  information.  For  example,  in  the  case 
of  two-dimensional  linear  finite  element,  there  are  four  ordinary  integration  points,  as 
depicted  in  Figure  8.  Figure  8  presents  the  positions  of  four  integration  points  in  q-q 
normalized  coordinate  system.  Each  of  the  integration  points  represents  a  quarter  of  area 
of  the  element  as  shown  in  Figure  8.  When  the  element  subdivision  technique  is  used, 
the  q-q  normalized  coordinate  values  of  each  generated  integration  point  are  evaluated 
and  the  program  chooses  an  ordinary  integration  point  whose  strain  history  parameters 
are  used.  For  three-dimensional  case,  the  same  strategy  is  adopted. 


•  Material  model  A  or  B 

having  the  higher  priority  is  used. 

•  Strain  history  data  on  the  smaller 


Figure  7  A  region  of  multiply  assignmed  material  models  and  the  Gauss  points  in 

the  overlapping  elements 
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Figure  8  Gauss  points  in  an  element  (two-dimensional  quadrilateral  element) 
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2.  Damage  constitutive  model 


In  this  research,  we  assume  the  nucleation  and  growth  of  microvoids  in  polymeric 
material  and  the  dewetting  between  the  particles  and  polymeric  matrix  material.  Figure 
9  shows  the  schematic  view  of  damage  modes.  These  material  damages  reduce  the 
effective  area  of  material  section.  Thus,  the  material  stiffness  is  reduced  gradually  while 
progressive  material  damages  take  place. 

Two  kinds  of  damage  constitutive  laws  were  adopted.  They  are  described  in  this  chapter. 
The  isotropic  damage  model  following  Simo  and  Ju  [6,  7]  was  extended  to  separate 
dilatational/deviatoric  damage  model  that  can  account  for  the  dilatational  and  deviatoric 
components  of  damages  separately.  The  separate  dilatational/deviatoric  damage  model 
is  especially  powerful  when  the  deformation  of  the  material  under  an  ambient  pressure 
is  analyzed. 

The  dilatational  part  of  material  damage  is  mainly  due  to  the  nucleation  and  growth  of 
microvoids.  It  is  considered  that  under  negative  hydrostatic  stress,  the  nucleation  and 
growth  of  microvoids  do  not  occur.  The  dilatational  damage  model  is  assumed  to  be 
induced  when  the  hydrostatic  stress  is  positive. 


I  t  t  t 


Matrix  damage  (microvoids,  etc.) 

iiii 

Figure  9  A  schematic  view  of  two  kinds  of  damage  modes;  the  growth  of 
microvoids  in  matrix  material  and  dewetting  between  the  particles  and  matrix 
material 

2.1  Isotropic  damage  constitutive  model 
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First,  we  describe  the  isotropic  damage  theory  by  following  Simo  and  Ju  [6],  In  Simo 
and  Ju  [6],  the  effective  stress  concept  and  the  hypothesis  of  strain  equivalence  are 
described.  When  they  are  applied  to  the  case  of  elastic  damage,  the  elastic  potential 

energy  of  damaged  material,  in  terms  of  the  strains  s and  the  damage  parameter  d 
is  written  to  be: 


y/(£U,d)  =  {\-d)f/0(£U) 


(12) 


where  y°{£u)  is  the  elastic-potential  function  for  virgin  material  that  is  written  to  be: 


¥  °  (£kt )  -  ^  Cijki  £ij  £  k( 


(13) 


Just  like  equivalent  stress  concept  in  the  theory  of  plasticity,  a  scalar  parameter  r  that 
measures  the  magnitude  of  stress/deformation  is  introduced. 

f  =  ^2i/’(sk,)  (14) 

Criteria  for  the  damage  evolution  are  written  as  follows.  First,  the  effective  stress  needs 
to  reach  the  critical  value. 

g{* {£U  ) -  r(d ))  =  {£k(  )-r{d)  =  0  (15) 

where  r(d )  is  the  function  of  the  damage  parameter  d  .  The  relationship  between  r(d) 
and  d  needs  to  be  determined  based  an  experimental  data.  The  material  damage 
progresses  when  the  damage  parameter  increases.  Therefore,  when  the  damage  is 
ongoing,  the  time  derivative  of  the  damage  parameter  d  must  be  positive. 

d>  0  (16) 

d  is  determined  through  the  evolution  law,  as: 

d  =  rH{r  ,d)=  t(£]c()h(t  ,d)  (17) 
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The  constant  H(r,d)  characterizes  the  progress  of  damage  with  respect  to  the  effective 
stress  like  term.  Equations  (15)  and  (16)  must  be  satisfied  when  the  progressive  damage 
takes  place.  It  is  noted  that  the  constant  H(f,d)  is  always  positive.  Therefore,  the 
equivalent  stress-like  term  must  increase  to  satisfy  equation  (17),  while  the  material 
damage  is  in  progress. 

Since  the  elastic  potential  energy  is  assumed  in  the  form  of  equation  (1),  the  stresses 
a ij  are  written  in  the  following  form. 


G ij  ~  (l  “  d  )Cijki8kf  (18) 

The  rates  of  stresses  can  be  expressed  by  differentiating  both  the  sides  of  equation  (12), 
as: 

Gij  —  (l  —  d)Cijki£k(  —  dCijk(£k(  —  (l  —  d)Cjjkf£  —  d&jj  (19) 


where  a ij(=Cijk(£k()  are  the  stresses  when  the  virgin  material  is  assumed  for  the  same 
strains.  From  equation  (17),  the  rate  d  of  damage  parameter  d  can  be  derived  to  be: 


d  =  tH(r,d)=  h(t  ,d)^-^-£ke 

T 


Hi^d) 


G kt 8 kt 


(20) 


H(z,d)  characterizes  the  rate  of  damage  parameter  d  with  respect  to  the  rate  of  the 
effective  stress-like  parameter  t . 

Substituting  equation  (20)  in  equation  (19),  we  arrive  at: 


-  (l  -  d]Cijke£ke  -  dcr°j  -  (l  -d)CijM£kt  ~ 


H(r,d) 


GijGUskl 


~  D!jk('8kf 


(21) 


where 
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Dijk(  ~  (l  -d)Cijki  - 


(22) 


H(r,d )  _0_0 

-  <JijaM 

D'jt,  are  the  tangent  moduli  for  the  isotropic  damage  material  and  have  the  major  and 
minor  symmetries,  i.e.  ofk,  =  ,  d$(  =  D$k  and  D,™  =  Dfkt . 


2.2  Separate  Isotropic/deviatoric  damage  model 

The  damage  evolution  of  some  class  of  materials  are  more  sensitive  to  hydrostatic 
pressure  stress  than  shear  stresses.  A  typical  scenario  is  in  a  material  containing 
microvoids.  Microvoids  grow  under  applied  positive-hydrostatic  pressure  stress.  But 
they  do  not  grow  under  negative-hydrostatic  pressure  stress.  The  growths  of  the  voids 
are  assumed  to  be  less  sensitive  to  the  shear  (deviatoric)  stresses  than 
positive-hydrostatic  stress. 

One  way  to  model  such  material  is  to  separate  the  contributions  of  isotropic  and 
deviatoric  stresses  to  the  damage  growth.  To  do  so  with  a  very  simple  model,  we 
separate  the  damage  parameter  d  into  the  dilatational  (volumetric)  and  deviatoric  parts. 
Hence,  the  elastic  potential  energy  function  [equation  (12)]  is  modified  and  is  written  to 
be: 


w{£id  ,dv 'd  o)  ~  -  dy)y/y  {ekk  )+{\-d  D  )//  °D 


(23) 


where  ekk  and  s\j  are  the  volumetric  and  deviatoric  strains.  The  functions  Vv  and 
Y°d  are  defined  to  be: 


Yv 


=  ^K{skk)2  and 


YV  =  M£ij£ij 


(24) 


where  K  and  //  are  the  bulk  and  shear  moduli.  It  is  noted  here  that  under  a  constraint 
condition  dv  =  dD ,  the  separate  dilatational/deviatoric  damage  model  is  the  same  as  the 
isotropic  damage  model.  We  define  two  kinds  of  effective  stress-like  terms,  as: 
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rv  -  (su )  -  jK(Skk  )2  and  td  -  ^2 )  -  ^2 ~^e\j£\j 


(25) 


When  the  damages  are  in  progress,  the  following  conditions  must  be  satisfied. 

fv  =  rv  (dv )  and  akk  >  0  for  the  dilatational  damage  (26) 

fD=rv(dD )  for  the  deviatoric  damage  (27) 

The  evolution  equations  for  the  damage  parameters  are  written,  as: 

Kg 

dy  =  tvHv(tv  ,dy)  =  Hy(rV  ,dy)y[K£kk  =  Hy(rV  ,dy  )—====  £kk  (28) 

A HJ2 

do  =TdH  D{TD’di))  =  H  D{TD’dD)  ==  £ki  (29) 

^2/J£  pq£  pq 

The  constants  Hv(fv,dv)  and  HD(fD,dD )  govern  the  evolution  laws  for  the  damage 
parameters.  The  constants  are  the  functions  of  the  effective  stress-like  terms  fv  and 
fD  and  of  the  damage  parameters  dv  and  dD .  From  equation  (23),  we  can  write  the 
stresses,  as: 


dv{£kl’dy  ,dD) 


dsn 


=  (l-dv)SijKskk+2{l- 


dD  )m£'h 


(30) 


By  differentiating  both  the  sides  of  equation  (30)  and  making  use  of  equations  (28)  and 
(29),  we  can  write  the  rate  form  constitutive  equation,  as: 

d,y  =  (l  -  dv  ]SjjKskk  +  2(l  -  dD  )//£',)  -dvKskk  -  2dD/.i6'lJ 

=  (l-dv)SijK£kk  +  2(l  -  d D  )fi£]j  —  Sjj — -  —{Kskk)2su  —  (~L> ’  —  /J1  £'ij£'kl£'kl 

Tv  td 

=  DWD*kt 

(31) 

and, 
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(32) 


DwD  =  (1  -dv  )KSySkt  +  (1  -  dD  )M(sikSje  +  SjkSu ) 
Hv(zv,dv )/  \2oo  4HD(TD,dD)  2 

-  =  ,LL  £ij£kl 


Dyke°  are  the  tangent  moduli  that  relate  the  rate  of  stresses  to  those  of  strains.  It  is  noted 


that  DyUD  have  the  major  and  minor  symmetries,  i.e.  DykeD  =  DklijD  ,  D;jk,D  =  D\Jfk  and 


Dv-r>  =  dv-d 

u ijki  jikt 
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3.  Cohesive  zone  model  and  its  implementation  in  the  s-FEM  computer 
program 

The  cohesive  zone  model  was  implemented  in  the  s-FEM  computer  program.  In  present 
research,  the  cohesive  zone  model  accounts  for  separation  between  the  reinforcing 
particles  and  matrix  material.  In  this  section,  the  formulations  of  the  cohesive  zone 
model  (Chandra  [8],  Foulk,  Allen  and  Helms  [9],  etc.)  along  with  the  s-FEM  are 
presented. 

3.1  Cohesive  zone  model 

Interface  separation  law  is  characterized  by  the  cohesive  zone  model.  The  interface  in  a 
solid  is  depicted  in  Figure  10.  We  assume  that  such  interface  separation  takes  place 
while  progressive  dewetting  between  matrix  and  a  particle  occurs.  A  three-dimensional 
illustration  is  given  in  Figure  11.  Here,  we  define  two  kinds  of  coordinate  systems.  One 
is  the  global  coordinate  system  (x1,x2,x3)  which  is  fixed  in  the  space  and  the  other  is 
local  coordinate  system  (x1,x2,x3)  in  which  x,  and  x2  are  in  the  plane  of  the 
interface  and  x3  is  perpendicular  to  the  interface.  We  define  upper  and  lower  faces  as 
depicted  in  Figure  11.  The  positive  direction  of  the  x3  coordinate  and  the  normal 
direction  n  of  the  interface  are  defined  as  depicted  in  Figure  11.  Their  directions  are 
defined  from  the  lower  to  the  upper  surface.  Relative  displacements  A ufHZ  are  defined 
as  the  difference  of  the  displacements  ufHZ+  and  uf'HZ~  of  the  upper  and  the  lower 

faces,  respectively.  (•).  indicate  that  the  components  of  the  vector  are  written  in  the 
(x,  ,x2,x3)  local  coordinate  system. 

A  ufHZ  =ufHZ+ -ufHZ~  (33) 

By  using  the  relative  displacements  A ufHZ ,  we  define  a  dimensionless  parameter  k 
which  characterizes  the  separation  and  the  slips  of  the  interface,  as: 


k  =  . 


UufHZ 

St 


\2 


Uu$HZ 

a 


\2 


f  STi(3H7 
Sn 


(34) 


where  a  and  Sn  are  the  length  parameters  that  characterizes  the  distance  of  partial 
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separation  before  the  final  one.  Then,  an  energy  density  function  is  defined  in  terms  of 
X . 


wCHZ  =  WCHZ(x) 


(35) 


W  CHZ  =  WCHZ(/ l)  is  the  energy  potential  and  has  properties  of: 


WCHZ(; l  =  0)  =  0 
WCHZ(X^oo)=WCHZ 


(36) 


Separation 


w 


CHZ 


is  the  energy  required  to  open  unit  area  in  cohesive  zone. 


I  Separation 

Since  wCHZ  =  WCHZ(x )  represents  the  energy,  we  can  express: 

rsytjCHZ  aw  CHZ 

tttCHZ(i\  A  —CHZ  ,  OW  —CHZ  /  ,  0\ 

w  +^F)4"3  (a=u) 

Therefore,  cohesive  tractions  f,  are  written  to  be: 


(37) 


T  = 
1  a 


dW 


CHZ 


A&F) 

Thus,  we  can  express  7) ,  as: 


(a  =  1,2),  T3  = 


dW 


CHZ 


(38) 


F  = 


dwCHZ 

_  dWCHZ 

8X 

_dWCHZ 

A  u™Z 

1 

_  dWCHZ 

Au™z 

)  dX 

)  dX 

St 

'~x 

dX 

AS? 

dWCHZ 

_  dWCHZ 

dX 

_  dWCHZ 

Au3chz 

1 

_  dWCHZ 

AU3chz 

5(a«3otz) 

dX 

d(Au3CHZ) 

dX 

Sn 

~x 

dX 

*S? 

(39) 


(40) 


It  is  assumed  here  that  wCHZ  has  the  following  form: 

WCHZ(x)  =  Sna^f(x)  (41) 

Here,  (j„AX  is  the  maximum  normal  stress  and  the  maximum  value  of  f(x)  is  one. 
/(; l)  is  the  dimensionless  function  of  X  .  The  function  /(; t)  takes  its  values 
f(X  =  o)  =  0  and  f(x  -> «)  =  o .  By  using  equation  (41),  the  tractions  can  be  written,  as: 
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CHZ 


T<x  -I  .  —  CHZ 


8W 


dW 


CHZ 


dX 


dX 


XrfF)’ 


dWCHZ  1  .  _max  df(X)  4i7 

°n(Jn 


CHZ 

a 


dA 


St 


dX 


XSf 


f  o  A 
Sn 


*3  = 


st 

dW 


I  °n 
J 

CHZ 


MAX  Sit g 117  1  df(X) 


dW 


St  X  dX 

CHZ 


i^HZ\ 


dA 


-  CJ 


k-CHZ 
MAX  Au3 


d\Au^ 
1  df(X) 


dX 


dWCHZmCHZl  {i)  &-CHZ 

- ^ - —  0„l T„  - - 

o  a  n  n 


dX 


dX  Xdl 


X  dX 


(42) 


(43) 


In  many  literatures  such  as  [8]  and  [9],  equations  (42)  and  (43)  are  interpreted  to  be: 


Ta=aAAXF(X) 


A  u%HZ 


f,=a^F{x) 


Au3CHZ 


(44) 

(45) 


where, 


a  =  ■ 


(46) 


Or  we  can  write: 


<r  A -CHZ 

fa  = 

Ot  Ot 


=  CJ 


MAX 


f(x) 


aAHZ 


(47) 

(48) 


In  order  to  obtain  rate  form  tangent  moduli,  we  further  differentiate  equations  (47)  and 
(48)  with  respect  to  time. 


t  =^rX\8AF{X)Au™ 


1  d F(X) 


X  dX 


Sn  AuZHZAufHZ 


A  u 


A  CHZ 
1 


aAhzaAhz 


A  u 


A  CHZ 


Au™ZAuCHZ 


s,  sn 


A  -CHZ 

-Am3 


Likewise, 
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(50) 


T,.-L«r{Fm ,CHZ 


1  df(/l) 


X  dX 


—CHZ  a  —CHZ 


Au^  Au[ 


-Au 


-CHZ 


Au£HZAu£HZ 


Au 


-CHZ 


\TfCHZ\77CHZ 

A“3  A“3 — au?hz 


In  a  matrix  form,  we  can  write: 


\To\  =  C 7 


MAX 


f(/1 


—  00 


—  0 


0  0  — 


1  dF  (x) 

X  dX 


.LA 

st  st 

i 


AufHZ  AufHZ 


8n  atLhzatChz 


8t  St  §} 

1  Au3CHZAufHZ 


St 


SnSt 


1 


C  .-CHZ .-CHZ 
on  A Uy  Au 2 


St  8, 

1 


C  .-CHZ .-CHZ 

8n  Am2  Am2 


1  Au3CHZAu™Z 


St 


8n8, 


1  AufHZAu3HZ 


8t  8n 
-CHZ  a -CHZ 


st 


1  Au 


Au a 


1  Awo 


Au\ 


7- CHZ 


A^ 

A -CHZ 

Au2 

Au3HZ 


(51) 


Therefore,  the  tangent  matrix  in  the  rate  formulation  is  symmetric. 


Function  f(a)  have  some  choices  as  described  in  Chandra  [8].  In  present  analyses, 
following  Foulk,  Allen  and  Helms  [9],  we  let  f{x)  be: 


F{X)~( X-2X  +  X2)  (52) 

and  therefore,  AAAI  =  2  +  2X)  =  —(x  - 1) 

dX  4  v  ’  2  v  ’ 

In  order  to  show  the  significance  of  this  choice  of  function  f(x),  we  consider  the 
energy  per  unit  area  of  the  interface  due  to  an  opening  mode.  For  the  absence  of  AufHZ 
and  A u2hz  components,  we  can  compute  the  energy  Q,  as: 

Q  =  \5LtA^3HZ  )=  Jo  T38ndX  =  \l<jLAXSnX^  -  2X  +  X2]iX  =  ^<?nAX8n  (53) 


where 
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(54) 


f3  =  oTX  fU)Au^-  =  oTX  f(A)A  =  °nAX  -(i-2A  +  A2)a 
8 n  4 

and 

d(Au3HZ)=SndA  (55) 

From  the  expression  of  (54),  the  maximum  value  of  T3  is  o™AX  for  A  =  ^ . 

From  equations  (53),  (54)  and  (55),  one  can  state: 


•  o^AX  and  8n  determines  the  characteristics  in  normal  separation  (energy  and 
initial  stiffness  of  the  interface). 

aMAX 

•  — -  characterizes  the  initial  stiffness  of  the  interface. 


s 

In  addition,  by  examining  equation  (50),  one  can  find  that  —  characterizes  the 

8t 

stiffness  ratio  between  the  normal  and  the  tangential  separations  of  the  interface. 


Figure  10  Progressive  dewetting/interface  separation  between  a  particle  and  matrix 
material 
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Figure  11  The  cohesive  zone  and  associated  coordinate  systems  ((xl,x2,x3)  global 
and  (xvx2,x3)  local  coordinate  systems) 


3.2  Cohesive  zone  with  an  arbitrary  orientation 

In  the  previous  section,  the  cohesive  zone  model  was  described  in  the  (x, ,  x2 ,  x3 )  local 
coordinate  system.  In  this  section,  detailed  discussions  on  how  the  coordinate 
translations  from  the  local  to  the  global  coordinate  system  are  carried  out  are  given. 

We  assume  a  rate  form  formulation  and  the  rates  of  stresses  on  the  interface  are  written 
in  terms  of  the  rates  of  relative  displacements. 

Tj  =  kijAUj'HZ  (56) 


where  the  components  of  k,j 


are  found  from  equation  (51),  as: 
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(57) 


When  the  interface  is  oriented  in  an  arbitrary  direction,  matrix  k,j  in  the  (xl,x2,x3) 

local  coordinate  system  must  be  transformed  to  the  c  global  coordinate  system.  The 
directions  of  x,  and  x2  coordinate  axes  are  in  the  plane  of  the  interface.  x3  is 
perpendicular  to  the  interface. 

An  inter-element  face  constitutes  a  cohesive  zone  element.  The  local  normalized 
coordinates  (#-77),  the  local  coordinates  (xlfx2,x3)  and  the  global  coordinates 
(x1,x2,x3)  are  set  as  shown  in  Figure  12.  Unit  basis  vectors  with  respect  to  the 
(x\ ,x2,x3)  and  (x| , x2 , x3 )  coordinate  systems  are  designated  to  be  (ei,e2,e3)  and 
(el,e2,e3),  respectively. 

First,  we  let  the  directions  of  £  and  x,  coordinate  axes  be  the  same,  as  shown  in 
Figure  12.  Thus,  the  basis  vector  et  is  shown  to  be: 


(58) 


Then,  we  define  another  unit  vector  e2 ,  and  e2  and  e3  are  calculated  to  be: 


dxi  ) 
— ~ei 

/% , 

and  e?  =  \e\  x  e?  )/\ei  x  e ? 

/  drj 

U  \  A  L*  /J  |  X  L*  | 

(59) 
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The  components  ky  of  tangential  stiffness  of  the  cohesive  zone  are  transformed  to  the 
global  coordinate  system,  as: 


kij  =ei- [kkeeket )'ej  =ku(ek  'eifee  'ej) 


(60) 


ktj  show  the  stiffness  of  the  cohesive  zone  in  the  (xl,x2,x3)  global  coordinate  system. 


t  t  t  t  t  t 


\  \  \  \  \  \ 

Figure  12  The  (x, ,  x2 ,  x3 )  local  coordinate  system  and  its  unit  basis  vectors 
(eve 2,J3),and  (c -//)  normalized  coordinates 


3.3  S-FEM  formulation  with  the  cohesive  zone  model 

3.3.1  Formulations 

The  cohesive  zone  model  is  implemented  in  the  s-FEM  computer  program.  It  is 
assumed  that  the  cohesive  zones  are  only  implemented  in  the  local  finite  element 
models.  The  global  finite  element  model  is  not  responsible  for  the  cohesive  zone  model. 
In  the  s-FEM  formulation,  we  start  with  the  principle  of  virtual  work,  as  stated  below. 
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(61) 


Jn° 


dSui 

dxj 


Ein 


ijkl 


duk 

dxp 


dzr 


—  |qg  ^ibfdCl  +  \qqG  duitid{dat  j 

+  SSCHZ+  tfHZ+SufHZ+d(sCHZ+  )+  \SCHZ-  tfHZ-SlfHZ~d(sCHZ-  ) 


where  sCHZ+  and  SCHZ  are  both  the  faces  of  the  cohesive  zone.  tfHZ+  and  tfHZ 
are  the  rate  of  tractions  on  sCHZ+  and  sCHZ~ ,  respectively,  that  are  induced  by  the 
opening  of  the  cohesive  zone.  SitfHZ+  and  SufHZ~  are  the  variations  of  the 
displacement  rates  on  sCHZ+  and  sCHZ~ .  They  are  in  the  state  of  equilibrium  and 
therefore  we  can  write: 


:CHZ+  _  -CHZ- 
Ti 

Therefore,  we  have: 


(62) 


\SCHZ+  tfHZ+SufHZ+d(sCHZ+  )+  \SCHZ-  tfHZ-SufHZ-d(sCHZ-  ) 

= jscHZ+ifHz+(sufHz+  ~^Hz~y (schz+)  (63) 

The  traction  rates  ifHZ+  are  written  in  terms  of  the  rates  of  relative  displacements 
fcufHZ  on  the  cohesive  zone,  as  shown  in  equation  (56).  The  relative  displacements  are 
written  in  terms  of  jumps  in  the  displacements  across  the  cohesive  zone.  Thus,  tfHZ+ 
are  written  to  be: 


■CHZ+ 

li 


kjj  ( uj 


CHZ+ 


■ CHZ- I 

~UJ  ) 


(64) 


Where  ufHZ+  and  ufHZ  are  the  velocities  at  both  the  faces  of  the  cohesive  zone.  ktj 

are  the  tangent  stiffness  of  the  cohesive  zone,  as  described  in  previous  section.  The  right 
hand  side  of  equation  (64)  has  a  negative  sign  because  the  positive  direction  of  the 
surface  is  opposite  from  that  in  previous  section.  By  substituting  equation  (64)  and  by 
using  the  superposition  of  displacements  based  on  the  global  and  local  finite  element 
models,  we  can  write: 


dduf 
dx  i 


Eijkt 


duj 


dxo 


-dQG  + 


M  o&lj 

2.  — 

p= 1  dx j 


Eijkt 


du 


Lp 


dxf 


dnLp  =  Jn0  Sufbi  dQc  +  Jsn<3  A,f,d(oQG )  (65) 


29 


(66) 


r  8Su\p 


F  ^uk  AcyLP  |  f  8dui  F 

Eijkelh7dn  +lnLp^~Eijke^ 


duIrp 


k  dnLp+  z 


8Sukp  „  8ukq 

_ _ _  /7. „ _ 11 _ 


J£,2^  n  ^IJKI  n  -  1  J12^  n  ^llKl  ~  -  '  ^  Ji2^  “*  n  ^ijki  n  - 

OXj  OX/?  OXy  OX//  6y=l  OXj  OX/? 

+}5c«z+  (,sif^z+(^)  -  («^z+(^)  -  <si^z-(£^))d(5c/,z+(z^>)=  y„  Sukpbt  dnLp 


dnLp~Lq 


(p  =  i,2,3,---,m) 


Here,  we  assume  that  there  are  M  local  models  that  are  superposed  on  the  global  model. 
The  local  models  may  overlap  each  other  and  their  coupling  terms  appear  in  equation 
(66).  It  is  noted  that  though  the  cohesive  zones  of  different  local  model  regions  may 
intersect  each  other,  their  coupling  terms  do  not  appear  in  equation  (66).  That  is  because 
two  cohesive  zones  make  a  line  of  intersection  when  they  intersect  and  the  line  has  a 
zero  area. 

3.3.2  Cohesive  zone  element 

When  the  cohesive  zone  is  implemented  in  the  s-FEM  computer  program,  we  introduce 
a  concept  of  cohesive  zone  element.  The  cohesive  zone  element  behaves  like  an 
interface  element.  Figure  13  shows  the  cohesive  zone  element.  For  an  illustrative 
purpose,  the  cohesive  zone  element  in  Figure  13  a  small  thickness.  In  actual  model,  the 
thickness  is  zero. 


The  terms  that  are  related  with  the  cohesive  zone  element  in  equation  (66)  as  expanded 
further,  as: 


\schz+  (sufHZ+iLp)  - &fHZ-{Lp%  (MfZ+(^  -  ^CHZ-iLp^CHZ^ 

=  \schz+  8u^Lp\ufz^  d(sCHZ+^)-  jscHZ+  SuCHZ^Lp\u<rHZ-(Lph (sCHZ<Lp))  (67) 

-  jscHz+  SufHZ~(Lphy<jHZ+(Lp^d{sCHZ+(Lp',)+  Ischz+  SufHZ~^Lphy^HZ~^d(sCHZ+^) 


Discretization  procedures  for  the  first  term  in  the  right  hand  side  of  equation  (67)  are 
described,  as  an  example.  Element  stiffness  matrix  for  a  cohesive  zone  element 
schz+(lp )  js  derjveci  as  follows.  The  displacements  within  sCHZ+(Lp^  are  expressed  by 
using  the  shape  functions  N1  (i  =  1,2, 3, 4)  where  the  CHZ  element  is  assumed  to  be 
linear  quadrilateral  element. 

«,-=M{«(l~4))  (68) 
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The  detail  of  equation  (68)  can  be  written  to  be: 


Ui 

N1 

0 

0 

N2 

0 

0 

N3 

0 

0 

N4 

0 

1 - 

o 

u2 

>  = 

0 

Nl 

0 

0 

N2 

0 

0 

N3 

0 

0 

N4 

0 

u3 

0 

0 

Nl 

0 

0 

N2 

0 

0 

N3 

0 

0 

N4 

u{ 

u2 

u\ 

u2 

u\ 

•2 

u3 

u\ 

«2 

il\ 

•4 

Ui 

•4 

u2 

■4 

u3 


(69) 


The  superscripts  in  equation  (69)  designate  the  nodal  numbers  (1,  2,  3  or  4). 

The  variations  of  the  velocities  can  also  be  written  in  the  same  manner.  The  stiffness  of 
the  CHZ  can  be  written  in  a  matrix  from,  as: 


:ll 

kl2 

kl3 

21 

k22 

k23 

31 

k32 

k33 

Thus,  we  have: 


(70) 


^CHZ+(LP)kij.CHZ+(LP)d^CHz4Lp^=  \^CHZ+(lP)^~4 


{aMf  jsc»z+^F[fcI^v]d(5c//z+^))Ji(i~4)} 


(71) 


The  other  terms  are  expressed  in  the  same  manner,  as: 


JSC*Z+  &fffZ+^)^Ji5’ffZ-^)d(5CTZ+^))=  {&M)}r  \kchz+(lp)^) 


{s(1~4)f  j$CHZ+  [/Vf  F][A']dkC//Z+^)fr(5-8)} 


(72) 
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=  jsi(5~8)|r 


(73) 


Jsc«z+  ^f//Z-^)^li5://Z+^>d(^e^-^))=  (si(5~8>f  [^C7^Z+(^)](i(5-8) 


:{»(5-8)fjsa«z+[^nifclAr]d(sCflZ^)^5-8)} 


(74) 


Thus,  the  term  of  the  CHZ  element  can  be  written  as  its  final  form,  as: 


jsCHZ+  [sufHZ+{^  - &CffZ-(LP)^(.CHZ+(LP)  _ &CHZALp)y(sCHZ+(LP)} 


=  {Si(1~4)[[^^Z+(L/))|i(l~4)}_  ^~4)J[KCHZ+(LP)p5~8) 

-  {si(5'8)^ [^CHZ+(LP)|i(l~4)}+  {S,0~4)jr^CWZ+(rP)|/(l--4)j 

KCHZ+(Lp)  _  KCHZ+(Lp ) 

_  KCHZ+(Lp)  KCHZ+(Lp ) 


^~8)  f 


(75) 


where 


^CffZ+(Lp)j=  |^wz+  [A?F[/f][,v]cl(1S’CWZ+(4')) 


(76) 


t  t  t  t  t  t 


nun 


Figure  13  A  schematic  view  of  a  cohesive  zone  element 
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4.  Material  damage  evolution  in  particulate  composite  materials: 
damage  constitutive  laws 

In  this  chapter,  we  discuss  about  the  mechanical  interactions  between  neighboring 
particles  and  the  distribution  and  the  evolution  of  material  damage.  We  assume  four 
kinds  of  materials;  (1)  the  isotropic  damage  material  and  the  separate 
dilatational/deviatoric  damage  material  with  the  constant  a  being  (2)  0.1,  (3)  0.99  and 
(4)  10.0.  For  the  separate  dilatational/deviatoric  damage  model,  the  distributions  of  the 
dilatational  and  deviatoric  damages  are  examined  separately. 

Matrix  material  is  assumed  to  undergo  the  damage.  The  sources  of  damage  evolutions 
are  due  to  the  nucleation  and  the  growth  of  microvoids.  The  stress-strain  curve  of  matrix 
material  is  postulated  to  be  that  of  Kwon  and  Liu  [10].  Although  the  uniaxial 
stress-strain  curve  of  Kwon  and  Liu  [10]  was  measured  for  a  polymeric  composite 
material,  it  is  used  to  model  matrix  material  in  this  investigation  since  there  is  no  other 
available  material  behavior. 

In  this  chapter,  we  discuss  about  (1)  how  the  relationship  between  the  effective 
stress-like  parameters  and  the  damage  parameters  are  obtained,  (2)  the  results  of 
analyses  for  unit  cell  models  are  described. 

4.1  Relationship  between  the  effective  stress-like  parameters  and  the  damage 
parameters. 

4.1.1  The  isotropic  damage  model 

The  damage  evolution  law  of  equation  (9)  has  a  scalar  function  H(f,d )  of  f  and  d  . 
H(f,d)  characterizes  the  damage  evolution  behavior  and  should  be  derived  from  a  set 
of  experimental  data.  In  this  section,  how  the  scalar  function  H(f,d)  is  determined  is 
described. 

First,  we  assume  that  we  have  a  stress-strain  curve  which  was  measured  in  an 
experiment.  Let  us  assume  that  we  have  a  uniaxial  stress-strain  curve,  as  shown  in 
Figure  14  Stress  is  expressed  by: 

a  =  {l-d)Es  (77) 
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Therefore,  the  damage  parameter  d  is  expressed  by  the  stress  and  strain,  as: 

d  =  1-—  (78) 

Es 

The  value  of  the  effective  stress-like  parameter  f  can  be  written  in  terms  of  the 
uniaxial  stress  and  strain,  as: 


Therefore,  once  the  uniaxial  stress-strain  curve  is  given,  relationship  between  the 
effective  stress-like  and  the  damage  parameter  can  be  obtained.  In  present  investigation, 
a  uniaxial  stress-strain  curve  is  approximated  by  a  series  of  piecewise  straight  lines,  as 
depicted  in  Figure  15.  In  Figure  15,  points  on  the  stress-strain  curve  ( cr1,s1 ),  (<t2 ,s2), 
(o-3,s3),  •••,  (a; ,£j),  •••  connect  straight  line  segments.  From  a  set  of 

data  (crj,^),  ((7 2 , £2 ) ,  {<t3,£3),  •••,  (o' j , £j ) ,  (<rM,ei+l) ,  •••  ,  we  can  compute 
(r, , ,  (f2,d2),  (i f3,d3 ),  •••,  (f, , dj ) ,  •••  by  using  equations  (78)  and 

(60).  Thus,  the  relationship  between  f  and  d  is  also  approximated  by  a  series  of 
piecewise  straight  lines.  The  function  H(f,d)  which  governs  the  damage  evolution  law 
of  equation  (20)  is  approximated  by: 

h(t, d)  =  4- «  — zr-  (dj  <  d  <  dI+i ,  r,-  <  r  <  r;+| )  (80) 

r  O+i  “0 

The  function  H(f,d )  whose  value  is  determined  by  equation  (80)  is  used  in  the 
stress-strain  relationship  of  equation  (21). 

We  then  extract  a  series  of  data  points  (o-, , ^ ) ,  (a2 , &-2 ) ,  (<j3 , £-3 ) , •  •  •  ,  (<t7,£7),  from 
the  experimental  stress-strain  curve  of  Kwon  and  Liu  [10],  as  shown  in  Figure  16.  As  an 
example,  the  relationship  between  the  uniaxial  strain  and  the  damage  parameter  for  the 
case  of  isotropic  damage  is  depicted  in  Figure  17.  Thus,  the  uniaxial  stress-strain  curve 
is  reconstructed  by  a  series  of  linear  approximations,  as  shown  in  Figure  18. 
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Slope:  E  (Young's  modulus) 


Fig.  14  Uniaxial  stress-strain  curve 


Slope:  E  (Young's  modulus) 


Fig.  15  Piecewise  linear  approximation  for  uniaxial  stress-strain  curve 


Y.W.  Kwon,  C.T.  Liu  /  Composites:  Part  B  32  (2001)  199-208 


Fig.  6.  Stress-strain  curve  of  a  uniform  composite  specimen. 


Figure  16  Uniaxial  stress-strain  curve  of  polymeric  particulate  composite 
material  that  was  given  in  Kwon  and  Liu  [10] 
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Figure  17  Relationship  between  the  strain  and  the  damage  parameter  for  the  isotropic 
damage  constitutive  model,  which  was  extracted  from  the  stress-strain  curve  of 
Figure  16. 


Figure  18  Uniaxial  stress-strain  curve  that  was  reconstructed  after  the  relationship 
between  the  effective  stress-like  and  the  damage  parameters  is  established  by  the 
proposed  procedures. 


4.1.2  The  separate  Isotropic/de viatoric  damage  model 

In  this  section,  the  damage  evolution  law  for  the  separate  isotropic/de  viatoric  damage 
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model  is  dealt  with.  As  seen  in  equations  (28)  and  (29),  there  are  two  different  damage 
evolution  laws  and  two  damage  variables  Hv  (fv ,  dv )  and  H D  (fD  ,dD).  In  order  to  fully 
determine  Hv  (fv , dv )  and  hd  (fD , dD ) ,  we  need  at  least  two  sets  of  stress-strain  curves. 
When  only  one  stress-strain  curve  is  available,  we  need  to  introduce  at  least  one 
additional  condition. 

We  assume  that  a  uniaxial  stress-strain  curve,  as  shown  in  Figure  15  or  16  is  available. 
As  an  additional  condition  to  uniquely  determine  the  damage  evolution  law,  we 
postulate  that  the  ratio  ( dv  jdD  )  between  the  dilatational  and  deviatoric  damage 
parameters  remains  to  be  the  same  during  the  uniaxial  deformation.  We  write: 

djy  —  ccdy  (81) 

Here,  a  is  a  positive  constant  ( 0  <  a  <  oo ).  When  a  equals  zero,  only  the  dilatational 
damage  is  present.  When  a  =  1 ,  the  magnitudes  of  dilatational  and  deviatoric  parts 
equal  each  other  and,  therefore,  this  case  is  very  similar  to  that  of  the  isotropic  damage. 
When  a  is  infinitely  large,  the  dilatational  damage  parameter  dv  is  zero.  Therefore, 
material  undergoes  deviatoric  damage  only. 

We  assume  uniaxial  deformation  in  x,  direction.  We  can  write  the  following 
statements. 


(Jjj  ^  0  (unknown),  <r22  =  cr33  =  <r12  =  cr23  =  cr31  =  0  (known) 
£|  |  i=-  0,  £'i 2  =  £'23  =  £'31  =  0  (known)  £^2  —  S'i'i  (unknown) 


(82) 


Therefore,  by  substituting  equations  (62)  and  (63)  in  equation  (82),  we  have: 


°il  -  {^-dy)K{s  11  +£22  +f33)+'^(l_a^v)//(2£li  -£22  “^33) 

2 

°22  =  (l  -  dv  )k{£\  1  +  £22  +  £33)+  —(l  -  adv)/j(2s22  -£33  -£n) 

2 

(Xu  =  (1_^vK(fll  +  e22  +£33 )  +  “(l_a4M2£33  ~£n  ~£22 ) 


(83) 


From  the  second  and  the  third  of  equation  (83)  and  s2 2  =  £33 ,  one  can  obtain: 
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(84) 


0  -  (l  -  dy  )K{t'\  |  +  £'22  +  f33  )+  “0  “  a^V  )/Ji£\  I  -£22) 


The  first  of  equation  (83)  and  equation  (84)  lead  to: 


-0 11 

£11  --Zu — ,  \~  +  £n 


2(l  -  adv  )/u 


(85) 


By  substituting  equation  (85)  in  the  first  of  equation  (83),  we  can  establish  a 
relationship  between  an  and  sx , ,  as: 


3(l  -  adv  Xl  -  dv  )/jK£i  |  -  j-i(l  -  adv )//  +  (l  -  dv  )K  JcT|  |  =0 
By  rearranging  equation  (86),  we  can  establish: 


(86) 


(xdw  + 


\aju  +  3^ 

^ll 

—  (l  +  CL  )  ftffy 

fi-  - ) 

1  9//A- 

*11 

J 

l  J 

=  0 


(87) 


Therefore,  for  a  give  set  of  variables  crn ,  sn  and  a  ,  we  can  solve  for  dv . 
Two  special  cases  ( a  =  0  and  a  =  1 )  are  presented  first.  When  a  =  0 ,  we  have: 


0.^+  °-//  +  3^^n 
1  9  juK  sn 


Wy  + 

fi-  ) 

1 

=  0 


and, 


d\/  —  1  —  - 


^11 
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When  a  =  1 ,  we  have: 
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(a  =  0) 


(91) 
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In  is  noted  that  there  are  two  solutions  that  satisfy  equation  (90)  and  they  are 

dv  =  1  -  and  dv  =  1 .  The  second  one  takes  a  constant  value  and,  therefore,  is  not  a 
Es  ii 

feasible  solution.  For  more  general  case  (a*  0,1 ),  we  solve  equation  (87)  for  dv ,  by 
letting: 

ccdy  +  Adv  +  B  =  0  (92) 

A=av± 3^n_(l  +  a)  and  B  =  l_*n_  (93) 

9  juK  sn  Es  ii 


Therefore,  the  value  dv  is  determined  to  be: 


dv  — 


-A+Va2  - 


4a# 


2a 


(94) 


In  equation  (94),  the  sign  associated  with  Va2  -4aB  needs  to  be  determined.  To 

determine  the  sign,  we  check  a  special  case  ( a  =  1 ).  By  letting  a  =  1 ,  in  equation  (94), 
we  have: 
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When  we  take  the  positive  sign,  equation  (95)  results  in: 
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This  is  a  constant  value  ( dv  =  1 )  and,  therefore,  is  inappropriate.  When  the  negative  sign 
is  assumed,  we  obtain: 
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This  result  is  the  same  as  the  case  of  isotropic  damage  material  (equation  (78)). 
Therefore,  we  choose  the  negative  sign  in  equation  (94). 


dv  —  ■ 


■A--\/a2  -4afi 
2  a 


Thus,  the  deviatoric  damage  parameter  dD  is  also  determined  to  be: 


dD  -  ady  =  a 


-A-Va2  -4 

2  a 


(98) 


(99) 


Once  uniaxial  stress-strain  curve  as  depicted  in  Figure  14  is  given,  we  correct  a  series 
of  points  (a2 , £'2 ) ,  (<t3,£-3),  •••,  (cr,- ,£-,),  (cr;+1 , s M ) ,  •••,  on  the  curve,  as 

shown  in  Figure  15.  Damage  parameters  ( dv\.  and  dD |. )  corresponding  to  the  stress 

and  strain  (cr, ,  )  are  obtained.  Thus,  by  using  equations  (82)  and  (83),  unknown 
strains  £-22  and  £-33  ( £-22  =  £-33 )  are  derived.  Therefore,  we  can  compute  the  effective 

stress-like  terms  tv  (=  4Kskk )  and  fD  (=  ) .  Thus,  we  have  a  series  of  data 

(dy\vdD\vfy\vfD\x)  ,  (4V|2,4D|2,FV|2,FD|2)  ,  (dv\ydD\3,TV\3,TD\3)  ,  •••  , 

(dv\.,dD\rTV\rTD\i),  (dv |.+1  ,dD | f+1 , FV |.+1  ,Fd|.+|),  ....  The  scalar  functions,  Hv (rv , dv ) 
and  HD(fD,dD )  are  derived  to  be : 
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5.  Material  damage  evolution  in  particulate  composite 
materials-cohesive  zone  model  (dewetting  between  the  particles  and 
matrix  material) 

5.1  Material  damage  evolution  in  matrix  material  under  the  influences  of 
mechanical  interactions  of  particles- 1  (Four  particle  model  without  ambient 
pressure) 

In  this  section,  we  try  to  present  how  the  damage  zone  develop  when  particulate 
composite  materials  are  loaded.  First,  we  consider  a  simple  problem  such  that  a  block 
containing  four  particles  is  subject  to  tension,  with  or  without  ambient  pressure. 

In  Figure  19,  a  model  that  contains  four  particles  is  presented.  The  finite  element 
models  for  the  s-FEM  analysis  are  also  illustrated  in  Figure  19.  As  presented  in  previous 
chapter,  the  stress-strain  curve  of  Kwon  and  Liu  [10]  is  postulated  and  the  relationship 
between  the  effective  stress-like  parameters  and  the  damage  parameters  are  obtained. 
From  the  results  of  this  simple  problem,  we  can  clarify  the  influences  of  particle 
arrangements  and  of  ambient  pressure  to  the  evolution  of  matrix  damage. 


Prescribed  Displacement 


Local  Model:  928  elements  1007  nodes 

Figure  19  The  problem  of  four  particles  with  and  without  the  ambient  pressure 


5.1.1  Isotropic  damage  material  without  ambient  pressure 
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The  distributions  of  stresses  and  damage  parameter  are  presented  in  Figure  20,  for  the 
levels  of  overall  strain  s  being  0.0632  and  0.1.  It  is  seen  that,  for  s  =  0.0632,  both 
tensile  and  transverse  stresses  concentrate  at  the  top  and  bottom  of  the  particles.  The 
magnitude  of  the  concentrations  seem  to  be  slightly  severer  in  the  gaps  between  the 
particles  than  the  other  locations.  Both  tensile  and  transverse  stresses  are  positive, 
resulsting  in  a  large  value  of  hydrostatic  stress.  Dmage  zones  start  to  develop  at  the  top 
and  bottom  of  the  particles.  In  the  gaps  between  the  particles,  the  damage  parameter  is 
slightly  larger  than  the  other  locations.  That  is  very  similar  to  the  case  of  stresses. 


For  s  =  0.1,  the  major  trends  are  similar  to  the  case  of  s  =0.0632.  However,  material 
damage  takes  place  almost  throughout  the  region  of  matrix  material. 
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Figure  20  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  isotropic  damage  without  ambient  pressure 


5.1.2  Separate  dilatational/deviatoric  damage  material  with  a  =  0.1  and  without 
ambient  pressure 
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In  the  case  of  the  separate  dilatational/deviatoric  damage  model,  the  dilatational  damage 
is  set  to  dominate  the  other.  In  Figure  21,  it  is  seen  that,  for  s  =  0.0632 ,  the  distributions 
of  the  tensile  and  transverse  stresses  are  similar  to  those  of  the  isotropic  damage  case. 
The  distributions  of  the  transverse  stress  look  differently.  That  is  because  the  ranges  are 
different.  The  dilatational  and  deviatoric  damage  zones  develop  at  the  top  and  bottom  of 
the  particles.  They  connect  the  particles  in  vertical  direction.  No  damage  develops  in 
horizontal  directions  from  the  particles.  As  expected,  the  dilatational  damage  parameter 
is  much  larger  than  the  deviatoric  one. 

When  the  level  of  overall  strain  s  increases  to  0.1,  the  damage  zones  enlarge  and  the 
most  of  part  of  matrix  material  suffer  from  the  material  damages.  The  dilatational 
damage  parameter  is  much  larger  than  the  deviatoric  one.  The  particles  constrain  the 
deformation  of  matrix  material  and  the  damage  variables  at  both  the  sides  of  the 
particles  are  almost  zero. 
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Figure  21  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  0.1  and  without  ambient  pressure 


5.1.3  Separate  dilatational/deviatoric  damage  material  with  a  =  0.99  and  without 
ambient  pressure 


In  this  case,  parameter  a  is  set  to  be  0.99.  The  distributions  of  tensile  and  transverse 
stresses  and  the  dilatational  and  deviatoric  damage  parameters  are  depicted  in  Figure  22. 
The  distributions  of  the  stresses  are  similar  to  afore  mentioned  cases.  Material  damages 
develop  at  the  top  and  bottom  of  the  particles  where  stresses  are  large.  As  shown  in 
Figures  22  (a-3),  (a-4),  (b-3)  and  (b-4),  the  levels  of  the  dilatational  and  deviatoric 
damage  parameters  are  similar.  Their  patterns  of  distributions  are  slightly  different.  The 
dilatational  damage  variable  at  the  sides  of  particles  is  almost  zero  and  such  regions 
connect  the  particles  in  horizontal  directions. 
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(a- 1)  Tensile  stress  (a-2)  Transverse  stress  (a-3)  Dilatational  (a-4)  Deviatoric 

damage  parameter  damage  parameter 

(a)  s  =  0.0632 


(b-1)  Tensile  stress  (b-2)  Transverse  stress  (b-3)  Dilatational  (b-4)  Deviatoric 

damage  parameter  damage  parameter 

(b)  £  =  0.1 

Figure  22  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  0.99  and  without  ambient  pressure 


5.1.4  Separate  dilatational/deviatoric  damage  material  with  a  =  10.0  and  without 
ambient  pressure 

In  this  case,  the  parameter  a  is  set  to  be  10.  It  is  assumed  that  the  deviatoric  damage 
mode  dominates  the  dilatational  one.  In  Figure  23,  the  distributions  of  the  stresses  and 
the  damage  parameters  are  shown. 

The  distributions  of  stresses  are  similar  to  afore  mentioned  cases.  The  level  of  the 
deviatoric  damage  parameter  is  much  larger  than  the  dilatational  one. 
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(a- 1)  Tensile  stress  (a-2)  Transverse  stress  (a-3)  Dilatational 

damage  parameter 

(a)  s  =  0.0632 


(a-4)  Deviatoric 
damage  parameter 


(b-1)  Tensile  stress  (b-2)  Transverse  stress  (b-3)  Dilatational  (b-4)  Deviatoric 

damage  parameter  damage  parameter 

(b)  £  =  0.1 

Figure  23  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  10.0  and  without  ambient  pressure 


5.1.5  Summary:  Separate  dilatational/deviatoric  damage  material  with  a  =  0.1, 
0.99  and  10.0,  and  without  ambient  pressure 

The  distributions  of  the  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  described  in  section  5. 1.1 -5.1.4.  The  stresses  are  large  at  the  top 
and  bottom  of  the  particles.  The  regions  of  high  stress  connect  the  particles  in  the 
vertical  direction.  At  the  sides  of  the  particles,  the  particles  constrain  matrix  material 
from  deformation.  Therefore,  the  damage  parameters  are  small  at  the  sides  of  the 
particles.  The  overall  trends  of  distributions  of  the  dilatational  and  deviatoric  damage 
parameters  are  very  similar.  Therefore,  it  is  summarized  that  the  choice  of  constitutive 
models  (the  isotropic  and  the  separate  dilatational/deviatoric  damage  models  with 
different  constant  a )  does  not  make  much  difference  in  the  results  when  no  ambient 
pressure  loading  is  applied  to  the  block. 
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5.2  Material  damage  evolution  in  matrix  material  under  the  influences  of 
mechanical  interactions  of  particles-2  (Four  particle  model  with  ambient  pressure) 

In  this  section,  the  same  four  particle  problems  are  solved  with  ambient  pressure  that  is 
applied  as  distributed  force  from  the  sides  of  the  block.  The  ambient  pressure  is  400kPa. 

Results  (the  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and 
deviatoric  damage  parameters)  are  visualized  for  overall  tensile  strain  s  being  0.1  and 
0.135.  Since  the  ambient  pressure  is  applied,  Poisson’s  ratio  effect  produces  negative 
tensile  stress.  Thus,  overall  tensile  strain  s  is  set  to  be  0.1  or  0.135  so  that  the  tensile 
stress  in  the  block  is  almost  the  same  as  the  cases  of  previous  section  (without  ambient 
pressure). 

5.2.1  Isotropic  damage  material  with  ambient  pressure 

In  Figure  24,  the  distributions  of  tensile  and  transverse  stresses  and  the  isotropic 
damage  parameter  are  depicted.  Major  trends  in  the  distributions  of  tensile  stress  are 
similar  to  those  without  the  ambient  pressure.  Also,  major  trends  in  the  distributions  of 
transverse  stress  are  also  analogous  to  those  without  the  ambient  pressure,  except  for 
that  the  value  of  stress  is  about  -400  kPa  at  distant  points  from  the  particles  whereas  it  is 
about  zero  for  the  case  without  the  ambient  pressure. 

The  distributions  of  the  isotropic  damage  parameter  are  quite  different  from  those 
without  the  ambient  pressure.  When  the  block  is  subject  to  tension  without  the  ambient 
pressure,  the  damage  zones  develop  at  the  top  and  bottom  of  the  particles  and  do  not 
appear  at  the  sides  of  the  particles.  However,  with  the  ambient  pressure,  the  damage 
zones  develop  at  the  sides  of  the  particles.  The  concentration  of  negative  stress  due  to 
ambient  pressure  occurs  at  the  sides  of  the  particles.  That  is  why  the  damage  zones 
develop  at  the  sides  of  the  particles. 

Development  of  damage  zone  due  to  the  growth  and  nucleation  of  microvoids  under 
negative  stress  are  not  likely  to  occur.  Therefore,  the  use  of  isotropic  damage 
constitutive  model  may  not  be  suitable  for  present  investigation. 
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Figure  24  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  isotropic  damage  with  ambient  pressure 


5.2.2  Separate  dilatational/deviatoric  damage  material  with  a  =  0.1  and  with 
ambient  pressure 

The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  depicted  in  Figure  25.  The  distributions  of  the  stresses  are 
similar  to  those  of  the  isotropic  damage  constitutive  law.  The  distributions  of  the 
damage  parameters  are  quite  interesting.  The  dilatational  damage  parameter  is  very 
small  compared  with  the  case  without  the  ambient  pressure.  The  distributions  of 
deviatoric  damage  parameter  are  almost  uniform  within  matrix  material.  The  level  of 
the  deviatoric  damage  parameter  is  small. 
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(a- 1)  Tensile  stress  (a-2)  Transverse  stress  (a-3)  Dilatational  (a-4)  Deviatoric 
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(b-1)  Tensile  stress  (b-2)  Transverse  stress  (b-3)  Dilatational  (b-4)  Deviatoric 

damage  parameter  damage  parameter 

(b)  s  =  0.135 


Figure  25  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  0.1  and  with  the  ambient  pressure 


5.2.3  Separate  dilatational/deviatoric  damage  material  with  a  =  0.99  and  with 
ambient  pressure 

The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  depicted  in  Figure  26.  It  is  seen  that  the  dilatational  damage 
parameter  is  much  smaller  than  the  deviatoric  one.  The  distributions  of  the  deviatoric 
damage  parameter  are  analogous  to  the  cases  of  the  isotropic  damage  constitutive  law 
with  the  ambient  pressure.  This  implies  that  the  development  of  damage  in  the  case  of 
the  isotropic  damage  law  is  governed  by  the  deviatoric  stresses. 
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Figure  26  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  0.99  and  with  the  ambient  pressure 


5.2.4  Separate  dilatational/deviatoric  damage  material  with  a  =  10.0  and  with 
ambient  pressure 

The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  depicted  in  Figure  27.  The  distributions  of  the  stresses  are 
similar  to  those  of  the  isotropic  and  the  separate  dilatational/deviatoric  damage  model 
with  the  other  values  of  parameter  a .  As  the  parameter  a  is  set  to  be  10.0,  the 
dilatational  damage  parameter  is  much  smaller  than  the  deviatoric  one.  The  distributions 
of  the  deviatoric  damage  parameter  look  similar  to  previous  cases  that  are  of  the 
isotropic  and  the  dilatational/deviatoric  damage  with  a  =  0.1  and  a  =  0.99. 
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Figure  27  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  10.0  and  with  the  ambient  pressure 

5.2.5  Summary:  Separate  dilatational/deviatoric  damage  material  with  a  =  0.1, 
0.99  and  10.0,  and  with  the  ambient  pressure 


Results  with  the  ambient  pressure  revealed  that  when  the  isotropic  damage  model  is 
used,  the  progressive  material  damage  occurs  even  under  negative  hydrostatic  stress. 
This  contradicts  with  our  basic  postulation  that  the  damage  in  matrix  material  is  due  to 
the  growth  and  nucleation  of  micro  voids.  These  should  not  be  sensitive  to  the  deviatoric 
stresses.  Therefore,  it  is  plausible  to  use  the  separate  dilatational/deviatoric  damage 
model  with  a  small  value  of  constant  a  . 


From  present  sets  of  analyses  we  are  unable  to  suggest  the  most  suitable  value  for  a  . 

5.3  Material  damage  evolution  in  matrix  material  under  the  influences  of 
mechanical  interactions  of  particles-3  (Nine  randomly  distributed  particle  model 
with  and  without  ambient  pressure) 


51 


Following  the  four  particle  problems,  we  solved  a  problem  with  nine  randomly 
distributed  particles.  The  problem  configuration  is  depicted  in  Figure  28.  Again,  we 
adopt  the  isotropic  damage  model  and  the  separate  dilatational/deviatoric  damage 
constitutive  model  with  «■  =  0.1,  a  =  0.99  and  a  =  10.0 .  Analyses  are  carried  out  with 
and  without  the  ambient  pressure. 


Local  Model:  928  elements  1007  nodes 

Figure  28  The  problem  of  none  randomly  distributed  particles  with  and  without  the 

ambient  pressure 

5.3.1  Isotropic  damage  material  without  ambient  pressure  (Nine  randomly 
distributed  particle  model) 

The  distributions  of  tensile  and  transverse  stresses  and  the  damage  parameter  are 
presented  in  Figure  29.  Trends  are  similar  to  the  case  of  four  regularly  distributed 
particles. 

It  is  seen  in  Figures  29  (a-3)  and  (b-3)  that  the  zone  of  large  value  of  the  damage 
parameter  connects  the  particles  when  they  are  close  to  each  other.  The  stresses 
concentrate  at  the  top  and  bottom  of  the  particles.  The  transverse  stress  at  the  sides  of 
the  particles  is  negative  and  the  zones  of  negative  stress  connect  the  particles  when  the 
distances  between  particles  are  small. 
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(b)  J  =  0.1 


Figure  29  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  isotropic  damage  without  ambient  pressure  (Nine  randomly 

distributed  particle  model) 


5.3.2  Separate  dilatational/deviatoric  damage  material  with  a  =  0.1  and  without 
ambient  pressure  (Nine  randomly  distributed  particle  model) 


The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  presented  in  Figure  30.  The  overall  trends  are  similar  to  that  of 
four  regularly  distributed  particle  case.  Again,  the  concentration  of  the  transverse  stress 
is  not  very  clear. 


The  value  of  the  dilatational  damage  parameter  is  much  larger  than  that  of  the  deviatoric 
one. 
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(a- 1)  Tensile  stress  (a-2)  Transverse  stress  (a-3)  Dilatational  (a-4)  Deviatoric 

damage  parameter  damage  parameter 
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(b-1)  Tensile  stress  (b-2)  Transverse  stress  (b-3)  Dilatational  (b-4)  Deviatoric 

damage  parameter  damage  parameter 

(b)  £=0.1 

Figure  30  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  0.1  and  without  the  ambient  pressure  (Nine  randomly  distributed  particle  model) 


5.3.3  Separate  dilatational/deviatoric  damage  material  with  a  =  0.99  and  without 
ambient  pressure  (Nine  randomly  distributed  particle  model) 

The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  depicted  in  Figure  31.  The  concentrations  of  the  transverse 
stress  at  the  top  and  bottom  of  the  particles  appear  more  clearly  than  the  case  of  a  =  0.1 . 
The  tensile  stress  concentrates  at  the  top  and  bottom  of  the  particles  as  seen  in  all  the 
other  cases.  The  values  of  the  dilatational  and  deviatoric  damage  parameters  are  similar. 
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(a- 1)  Tensile  stress  (a-2)  Transverse  stress  (a-3)  Dilatational  (a-4)  Deviatoric 

damage  parameter  damage  parameter 

(a)  s  =  0.0632 


(b-1)  Tensile  stress  (b-2)  Transverse  stress  (b-3)  Dilatational  (b-4)  Deviatoric 

damage  parameter  damage  parameter 

(b)  £  =  0.1 

Figure  3 1  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  0.99  and  without  the  ambient  pressure  (Nine  randomly  distributed  particle  model) 


5.3.4  Separate  dilatational/deviatoric  damage  material  with  a  =  10.0  and  without 
ambient  pressure  (Nine  randomly  distributed  particle  model) 

The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  depicted  in  Figure  32.  As  seen  in  Figures  32  (a-2)  and  (b-2),  the 
concentrations  of  the  transverse  stress  at  the  top  and  bottom  of  the  particles  are  clearly 
seen.  The  value  of  the  dilatational  damage  parameter  is  much  smaller  than  that  of  the 
deviatoric  one. 
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Figure  32  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  10.0  and  without  the  ambient  pressure  (Nine  randomly  distributed  particle  model) 


5.3.5  Isotropic  damage  material  with  ambient  pressure  (Nine  randomly  distributed 
particle  model) 

The  distributions  of  tensile  and  transverse  stresses  and  the  damage  parameter  are 
depicted  in  Figure  33.  As  seen  in  the  case  of  regularly  placed  four  particle  problem, 
progressive  material  damage  at  the  sides  of  the  particles  are  seen.  They  are  due  to 
negative  transverse  stress. 
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Figure  33  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  isotropic  damage  with  the  ambient  pressure  (Nine  randomly 

distributed  particle  model) 


5.3.6  Separate  dilatational/deviatoric  damage  material  with  a  =  0.1  and  with  the 
ambient  pressure  (Nine  randomly  distributed  particle  model) 


The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  presented  in  Figure  34.  The  distributions  of  the  dilatational  and 
deviatoric  damage  parameters  are  quite  different  from  that  of  the  damage  parameter  of 
the  isotropic  damage  model.  Also,  they  look  differently  from  the  case  without  the 
ambient  pressure.  It  seems  that  the  ambient  pressure  weakens  the  magnitude  of 
hydrostatic  stress  and  prevents  the  dilatational  damage  zone  from  enlarging  in  matrix 
material.  The  deviatoric  damage  parameter  looks  more  uniform  than  the  case  without 
the  ambient  pressure. 
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(a- 1)  Tensile  stress  (a-2)  Transverse  stress  (a-3)  Dilatational  (a-4)  Deviatoric 

damage  parameter  damage  parameter 

(a)  J  =  0.1 


(b-1)  Tensile  stress  (b-2)  Transverse  stress  (b-3)  Dilatational  (b-4)  Deviatoric 

damage  parameter  damage  parameter 

(b)  s  =  0.135 

Figure  34  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  0.1  and  with  the  ambient  pressure  (Nine  randomly  distributed  particle  model) 


5.3.7  Separate  dilatational/deviatoric  damage  material  with  a  =  0.99  and  with  the 
ambient  pressure  (Nine  randomly  distributed  particle  model) 


The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  presented  in  Figure  35.  The  levels  of  the  dilatational  damage 
parameter  are  lower  than  the  case  of  a  =  0.1.  In  the  distributions  of  the  deviatoric 
damage  parameter,  the  concentrations  of  the  damage  parameter  are  seen  not  only  at  the 
top  and  bottom  of  but  also  at  both  the  sides  of  the  particles.  It  seems  that  the 
compressive  stress  due  to  the  ambient  pressure  induces  the  evolution  of  the  deviatoric 
damage. 
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(a- 1)  Tensile  stress  (a-2)  Transverse  stress  (a-3)  Dilatational  (a-4)  Deviatoric 

damage  parameter  damage  parameter 
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(b-1)  Tensile  stress  (b-2)  Transverse  stress  (b-3)  Dilatational  (b-4)  Deviatoric 

damage  parameter  damage  parameter 

(b)  s  =  0.135 

Figure  35  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  0.99  and  with  the  ambient  pressure  (Nine  randomly  distributed  particle  model) 


5.3.8  Separate  dilatational/deviatoric  damage  material  with  a  =  10.0  and  with  the 
ambient  pressure  (Nine  randomly  distributed  particle  model) 


The  distributions  of  tensile  and  transverse  stresses  and  the  dilatational  and  deviatoric 
damage  parameters  are  presented  in  Figure  36.  The  major  trends  are  the  same  as  the 
case  of  a  =  0.99 .  The  value  of  the  deviatoric  damage  parameter  is  slightly  larger  than 
that  of  a  =  0.99 .  The  value  of  the  dilatational  damage  parameter  is  much  smaller  than 
those  of  a  =  0.10  and  a  =  0.99  . 
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(a- 1)  Tensile  stress  (a-2)  Transverse  stress  (a-3)  Dilatational  (a-4)  Deviatoric 

damage  parameter  damage  parameter 

(a)  £=0.1 


(b-1)  Tensile  stress  (b-2)  Transverse  stress  (b-3)  Dilatational  (b-4)  Deviatoric 

damage  parameter  damage  parameter 

(b)  s  =  0.135 

Figure  36  The  distributions  of  the  tensile  and  transverse  stresses  and  the  damage 
parameter  for  the  case  of  the  separate  dilatational/deviatoric  damage  model  with 
a  =  10.0  and  with  the  ambient  pressure  (Nine  randomly  distributed  particle  model) 

5.3.9  Summary:  Separate  dilatational/deviatoric  damage  material  with  a  =  0.1, 
0.99  and  10.0,  and  with  and  without  the  ambient  pressure  (Nine  randomly 
distributed  particle  model) 

From  the  sets  of  analyses  presented  in  section  5.3,  it  is  clearly  seen  that  the  major  trends 
in  the  case  of  randomly  distributed  particles  are  the  same  as  those  in  the  four  particle 
problem.  However,  it  is  seen  that  the  damage  zones  connect  the  particles  depending  on 
their  distances.  The  most  important  finding  in  section  5.3  is  as  follows.  The  results  seem 
to  be  somewhat  reasonable  with  all  the  constitutive  models,  i.e.,  the  isotropic  and  the 
separate  dilatational/deviatoric  damage  model  with  a  =  0.1,  a  =  0.99  and  a  =  10.0 
when  the  ambient  pressure  is  not  applied.  However,  when  the  ambient  pressure  is 
applied  to  the  block,  the  deviatoric  damage  zones  develop  at  both  the  sides  of  the 
particles,  where  the  hydrostatic  stress  is  negative.  In  present  study,  we  assume  that 
material  damage  is  due  to  the  growth  and  nucleation  of  microvoids  in  matrix  material. 
Such  material  damage  is  not  considered  to  take  place  under  negative  hydrostatic  stress. 
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Therefore,  it  is  appropriate  to  use  the  separate  dilatational/deviatoric  damage 
constitutive  model  with  a  small  number  of  constant  a  .  When  the  constant  a  is  small, 
the  contribution  of  the  deviatoric  damage  is  small  relative  to  that  of  the  dilatational 
damage. 

5.4  Material  damage  evolution-dewetting  between  the  particles  and  matrix 
material  (cohesive  zone  model) 

In  this  section,  some  preliminary  results  are  presented.  Due  to  the  softening  nature  of 
the  cohesive  zone  model,  s-FEM  analyses  tend  to  become  unstable.  Thus,  we  could  not 
use  iterative  solver  to  solve  for  the  system  of  linear  simultaneous  equations.  Thus,  we 
used  a  direct  solver  (skyline  solver)  in  the  analyses  of  section  5.4,  although  it  is  not  very 
efficient  in  terms  of  its  memory  usage.  It  is  noted  that  for  all  the  other  analyses,  an 
iterative  equation  solver  (Conjugate  Gradient  Method)  is  used.  Because  of  this  problem, 
we  could  place  as  many  as  four  particles. 

Since  the  constants  8n ,  S,  and  <j„ax  in  the  cohesive  zone  model  are  not  known, 
postulated  values  are  used  in  the  analyses.  We  let  8n  and  8,  be  0.01  which  is  much 
smaller  than  the  size  of  the  particle  in  out  problem.  o™AX  is  varied  and  from  small  to 
large  values  relative  to  the  nonlinear  stress-strain  curve  of  Figure  16,  i.e.,  <yAMX  =250 
kPa,  500  kPa,  1000  kPa  and  2000  kPa. 

5.4.1  One  particle  model,  without  matrix  damage  (without  the  ambient  pressure) 

In  this  section,  the  results  of  one  particle  problem  that  is  shown  in  Figure  37  are 
presented.  The  particle  is  placed  in  the  block  made  of  elastic  material.  The  constants  for 
the  cohesive  zone  model  are  postulated  to  be: 

Sn  =  0.01 ,  st  =  0.01 ,  a^AX  =250  kPa,  500  kPa,  1000  kPa  or  2000  kPa  (101) 

In  Figure  38,  the  tensile  stress  distributions  in  the  central  section  of  the  block  are  shown 
for  the  case  of  <j„ax  =250  kPa.  In  Figure  38  (a),  the  stress  distribution  shows  stress 
concentrations  at  the  top  and  bottom  of  the  particle.  The  distribution  is  very  similar  to 
that  without  the  cohesive  zone  model.  This  means  that,  at  this  step,  the  cohesive  zone  is 
almost  rigid.  In  Figure  38  (b),  it  is  seen  that  the  stress  concentrations  at  the  top  and 
bottom  of  the  particle  are  weaker  than  those  in  Figure  37  (a).  The  cohesive  zone  opens 
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and  the  load  transmitting  capability  of  the  cohesive  zone  starts  decreasing. 

In  Figure  38,  the  distributions  of  the  tensile  stress  for  the  case  of  <j^ax  =500  are  shown. 
The  trends  are  similar  to  those  of  ajfAX  =250.  In  Figures  39  and  40,  the  distributions  are 
depicted  for  cr^AX  =  1 000  and  <j„ax= 2000,  respectively.  Except  for  Figure  39  (b),  the 
trends  are  very  similar  to  those  of  <j^ax  =250  and  a^AX  =500.  Due  to  the  softening 
nature  of  the  cohesive  zone  model,  the  analysis  suddenly  became  unstable.  That  is  why 
the  stress  distribution  Figure  39  (b)  looks  asymmetric. 


Prescribed  Displacement 


Local  Model:  928  elements  1007  nodes 

Figure  37  One  particle  problem  with  the  cohesive  zone  model 
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Figure  37  The  distributions  of  tensile  stress  for  <j^ax  =250  kPa 
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Figure  38  The  distributions  of  tensile  stress  for  a„AX  =500  kPa 
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Figure  39  The  distributions  of  tensile  stress  for  ajfAX  =1000  kPa 
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Figure  40  The  distributions  of  tensile  stress  for  <r^AX  =2000  kPa 

5.4.2  One  particle  model,  with  matrix  damage  (without  the  ambient  pressure) 


In  section  5.4.2,  in  addition  to  the  cohesive  zone  model,  we  assume  the  evolution  of 
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material  damage  in  matrix  material.  As  seen  in  preceding  sections,  the  separate 
dilatational/deviatoric  damage  model  with  a  small  value  of  constant  a  has  been 
recommended  to  appropriately  model  the  matrix  damage  due  to  the  growth  and 
nucleation  of  microvoids.  Thus,  in  this  section,  a  is  set  to  be  0.1.  The  stress-strain 
curve  of  Figure  16  is  assumed  to  model  the  behavior  of  matrix  material.  ajfAX ,  a 
constant  in  the  cohesive  zone  model  is  postulated  to  be  250  kPa,  500  kPa,  1000  kPa  or 
2000  kPa. 

In  Figure  41,  the  distributions  of  tensile  stress  and  of  the  dilatational  damage  parameter 
for  a„AX  =  250  kPa  are  presented.  It  is  considered  that  cr^AX  =  250  kPa  is  small 
compared  with  the  level  of  stress-strain  curve  of  Figure  16.  According  to  Figure  16,  the 
nonlinear  deformation  of  matrix  material  initiates  at  500  kPa.  The  value  a^AX  =  250 
kPa  is  smaller  than  the  stress  at  the  initiation  of  damage.  The  tensile  stress  distribution 
of  Figure  41  (b)  indicates  that  the  stress  concentration  weakens  due  to  the  opening  of 
the  cohesive  zone  at  the  top  and  bottom  of  the  particle.  Figure  41  (c)  shows  that  the 
dilatational  damage  parameter  is  very  small.  Therefore,  when  the  value  of  a™AX  is 
small  compared  with  that  of  damage  initiation,  the  opening  of  cohesive  zone  occurs 
before  the  material  damage  of  matrix  material  takes  place.  Therefore,  the  major  part  of 
the  material  damage  is  dewetting  between  the  particle  and  matrix  material. 

In  Figure  42,  the  distributions  of  tensile  stress  and  the  dilatational  damage  parameter  are 
shown  for  <r^AX  =  500  kPa.  In  this  case,  the  value  of  a]fAX  is  comparable  to  the  stress 
value  at  the  initiation  of  matrix  damage.  In  Figures  42  (a)  and  (c),  the  distributions  of 
tensile  stress  and  the  dilatational  damage  parameter  at  overall  strain  0.0441  are  depicted. 
At  this  stage,  the  damage  zone  has  not  developed  much.  The  tensile  stress  distribution  is 
similar  to  that  of  elastic  problem. 

In  Figures  42  (b)  and  (d),  the  distributions  of  tensile  stress  and  the  dilatational  damage 
parameter  are  presented  for  overall  strain  being  0.0632.  The  stress  distribution  indicates 
that  the  dewetting  between  the  particle  and  matrix  material  is  taking  place.  Figure  42  (d) 
shows  that  the  dilatational  damage  zones  are  developing  at  the  top  and  bottom  of  the 
particle.  In  this  case,  both  the  damage  modes  i.e.,  matrix  damage  and  dewetting  between 
the  particle  and  matrix  take  place. 

In  Figure  43,  the  distributions  of  tensile  stress  and  the  dilatational  damage  parameter  are 
presented  for  crjfAX  =  1000  kPa.  From  the  figures,  it  is  seen  that  the  cohesive  zone  does 
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not  influence  the  distributions  of  stress  and  the  damage  parameter.  The  cohesive  zone  is 
so  stiff  that  it  does  not  open  and  material  damage  due  to  the  dewetting  does  not  occur.  In 
Figure  44,  the  distributions  are  presented  for  <j„ax  =  2000  kPa.  They  are  almost 
identical  to  those  of  <j^ax  =  1000  kPa. 
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(a)  Tensile  stress  at  overall  strain=0.015  (b)  Tensile  stress  at  overall  strain=0.0373 
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(c)  Dilatational  damage  parameter  at  overall  strain=0.0373 


Figure  41  The  distributions  of  tensile  stress  and  of  the  dilatational  damage 
parameter  for  a^AX  =250  kPa. 
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Figure  42  The  distributions  of  tensile  stress  and  of  the  dilatational  damage 
parameter  for  <j„ax  =500  kPa. 
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(c)  Dilatational  damage  parameter  (d)  Dilatational  damage  parameter 

at  overall  strain=0.0632  at  overall  strain=0. 1 


Figure  43  The  distributions  of  tensile  stress  and  of  the  dilatational  damage 
parameter  for  a„AX  =1000  kPa. 
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Figure  44  The  distributions  of  tensile  stress  and  of  the  dilatational  damage 
parameter  for  a^AX  =2000  kPa. 


5.4.3  Four  particle  model,  with  matrix  damage  (without  the  ambient  pressure) 

In  this  section,  the  results  of  analyses  on  the  four  particle  model  are  presented.  The 
same  damage  constitutive  law  as  that  is  used  to  solve  the  one  particle  problems  is 
adopted  (the  separate  dilatational/deviatoric  model  with  a  =  0.1 ). 

In  Figure  45,  the  distributions  of  tensile  stress  and  the  dilatational  damage  parameter  for 
°nAX  =250  kPa  are  presented.  It  is  seen  in  the  figures  that  the  value  of  the  damage 
parameter  is  small  although  the  cohesive  zones  start  separating  and  that  the  stress 
concentration  at  the  top  and  bottom  of  the  particles  weakens  due  to  the  softening  of  the 
cohesive  zone,  i.e.,  dewetting. 

In  Figure  46,  the  distributions  of  tensile  stress  and  he  dilatational  damage  parameter  for 
anAX  =500  kPa  are  depicted.  In  this  case,  the  stress  at  the  maximum  stress  at  the 
cohesive  zone  and  that  at  the  initiation  of  matrix  damage  are  comparable.  It  is  seen  in 
Figure  46  (b)  that  the  stress  concentration  at  the  top  and  bottom  of  the  particles  weakens 
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due  to  both  the  damage  mechanisms.  The  damage  zones  connect  the  particles  aligned  in 
the  vertical  direction,  as  seen  in  Figures  46  (c)  and  (d). 

In  Figures  47  and  48,  the  distributions  of  tensile  stress  and  the  dilatational  damage 
parameter  are  presented  for  <j^ax  =1000  kPa  and  aAMX  =2000  kPa,  respectively.  The 
distributions  in  Figures  47  and  48  are  very  similar,  although  the  different  values  are 
given  to  the  parameter  a„AX .  In  both  the  cases  of  Figures  47  and  48,  the  cohesive 
zones  do  not  give  any  influences  to  the  deformations  of  the  composites.  The  cohesive 
zones  are  stiff  enough  that  they  do  not  open. 

As  summary,  we  can  conclude  that  depending  on  the  combination  of  stresses  at  the 
initiation  of  matrix  damage  and  at  the  separation  of  cohesive  zone.  In  present  analyses, 
the  stress-strain  curve  of  matrix  material  is  postulated  to  be  what  was  given  in  Kwan 
and  Liu  [10]  and  the  parameter  o-xlAX  of  cohesive  zone  model  that  governs  the  stress  at 
the  separation  of  cohesive  zone  is  varied.  When  <j„ax  is  set  to  be  250  kPa,  the  stress  at 
the  separation  of  cohesive  zone  is  smaller  than  that  at  the  initiation  of  matrix  damage.  In 
such  a  case,  interface  separation  at  the  cohesive  zone  is  the  major  damage  mechanism  of 
the  composite.  When  <j„ax  is  500  kPa,  the  stress  at  the  separation  of  cohesive  zone  is 
comparable  to  that  at  the  progress  of  matrix  damage.  Both  the  mechanisms  of  material 
damage  take  place.  When  aAMX  is  1000  kPa  or  2000  kPa,  the  interface  separation  at 
the  cohesive  zone  do  not  occur  at  all.  The  matrix  damage  is  the  major  damage 
mechanism. 
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(c)  Dilatational  damage  parameter  at  overall  strain=0.0411 
Figure  45  The  distributions  of  tensile  stress  and  of  the  dilatational  damage 
parameter  for  a„AX  =250  kPa  (four-particle  problem). 


L 


6. 52e+02 
5. 97e+02 
5. 43e+02 


4. 34e+02 
3. 79e+02 
3. 25e+02 
2. 70e+02 
j2. 16e+02 

I1.61e+02 

1.07e+02 


(a)  Tensile  stress  at  overall  strain=0.0411 


9. 78e+02 
9.01e+02 
8. 23e+02 


69e+02 

91e+02 

14e+02 


l 


59e+02 

82e+02 

04e+02 


(b)  Tensile  stress  at  overall  strain=0.0632 


(c)  Dilatational  damage  parameter 
at  overall  strain=0.04112 


(d)  Dilatational  damage  parameter 
at  overall  strain=0.0632 


Figure  46  The  distributions  of  tensile  stress  and  of  the  dilatational  damage 
parameter  for  a„AX  =500  kPa  (four-particle  problem) 
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Figure  47  The  distributions  of  tensile  stress  and  of  the  dilatational  damage 
parameter  for  o„AX  =1000  kPa  (four-particle  problem) 
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(c)  Dilatational  damage  parameter  (d)  Dilatational  damage  parameter 

at  overall  strain=0.0632  at  overall  strain=0.1 


Figure  47  The  distributions  of  tensile  stress  and  of  the  dilatational  damage 
parameter  for  o„AX  =2000  kPa  (four-particle  problem) 
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6.  Conclusions 


In  this  investigation,  the  deformation  and  damage  mechanisms  of  particulate  composite 
material  are  investigated.  To  do  so,  we  adopted  the  s-version  Finite  Element  Method 
(s-FEM).  By  using  the  s-FEM,  the  finite  element  models  for  the  structure  as  whole  and 
for  a  particle  and  its  vicinity  are  built  separately  and  are  superposed  each  other.  They  are 
called  the  global  and  the  local  finite  element  models.  Once  the  global  model  and  the 
local  model  are  built,  they  are  superposed  in  any  fashion  as  long  as  the  region  of  the 
local  model  is  included  in  that  of  the  global  model.  The  local  model  can  be  superposed 
on  the  global  model  repeatedly  and  they  can  also  overlap  each  other.  Thus,  as  seen  in 
this  report,  without  building  new  finite  element  models,  various  particle  arrangements 
can  be  modeled  without  creating  any  new  finite  element  models. 

In  this  research,  continuum  damage  constitutive  models  are  adopted.  They  are  based  on 
Simo  and  Ju  [6],  To  account  for  the  contributions  of  the  hydrostatic  and  deviatoric 
stresses  independently,  the  separate  dilatational/deviatoric  damage  constitutive  model 
was  adopted.  By  adjusting  a  constant  in  the  model,  we  can  model  various  combinations 
of  the  contributions  from  the  hydrostatic  and  deviatoric  stresses.  For  example,  when  the 
constant  a  is  set  to  zero,  purely  dilatational  damage  solid  is  modeled.  When  a  =  1 ,  the 
behavior  of  the  solid  is  similar  to  that  of  isotropic  damage  material,  a  ->  oo  results  in  a 
deviatoric  damage  material. 

Cohesive  zone  model  is  implemented  in  the  s-FEM  computer  program  to  account  for 
dewetting  between  the  particles  and  matrix  material.  Only  the  local  finite  element  model 
contains  the  cohesive  zone  as  the  global  model  does  not  account  for  the  details  of  the 
structure  such  as  the  particles.  The  cohesive  zone,  when  it  separates,  has  a  softening 
behavior.  The  iterative  equation  solver  (Conjugate  Gradient  Method)  that  is  adopted  in 
the  program  can  not  handle  the  softening  behavior.  Thus,  a  direct  solver  (Skyline 
method)  was  used  to  solve  the  problems.  The  direct  solver  requires  much  larger  memory 
space  than  the  iterative  one  does.  Therefore,  it  was  difficult  to  place  many  particles.  At 
most,  we  could  place  four  particles. 

The  results  presented  in  this  report  revealed  that  the  use  of  isotropic  damage  model  is 
not  appropriate  since  the  material  damage  may  progress  under  a  negative  hydrostatic 
stress.  The  use  of  the  separate  dilatational/deviatoric  damage  model  is  recommended 
with  the  hydrostatic  stress  having  the  major  contribution.  The  evolutions  of  progressive 
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damage  zones  are  presented.  When  the  particles  are  aligned  in  the  loading  direction,  the 
damage  zones  connect  the  particles.  Due  to  the  strong  constraint  from  the  particle, 
damage  zones  do  not  develop  at  the  sides  of  the  particle. 

Dewetting  between  the  particles  and  matrix  material  is  seen  to  initiate  at  the  top  and  the 
bottom  of  the  particles.  Dewetting  does  not  occure  at  the  sides  of  the  particles. 
Depending  on  the  strength  of  bond  between  the  particles  and  matrix  material,  dominant 
damage  mechanism  changes.  When  the  bond  is  very  strong,  i.e.,  the  value  of  constant 
crMax  is  large,  only  the  matrix  damage  occurs.  When  the  bond  is  very  weak,  i.e.,  crMax 
is  small,  dewetting  is  the  dominant  damage  mechanism.  When  the  strength  of  the  bond 
is  comparable  to  the  stress  at  the  initiation  of  matrix  damage,  both  the  damage 
mechanisms  take  place. 

Although  s-FEM  computer  code  has  been  developed  for  the  damage  analyses,  the 
material  constants  that  were  adopted  in  the  analyses  are  postulated  ones.  However,  some 
characteristic  behavior  of  particulate  composite  materials  with  progressive  damage  has 
been  revealed. 

In  summary,  during  present  research,  the  followings  have  been  accomplished. 

1)  The  separate  dilatational/deviatoric  damage  constitutive  model  was  proposed  along 
with  a  method  to  extract  its  material  parameters  from  one-dimensional  stress-strain 
curve.  The  damage  model  was  implemented  in  the  s-FEM  computer  program, 

2)  S-FEM  formulation  with  the  cohesive  zone  model  has  been  proposed  and  has  been 
implemented  in  the  s-FEM  program. 

3)  S-FEM  analyses  on  particulate  composite  materials,  with  assuming  the  matrix 
damage  and  dewetting  between  the  particles  and  matrix  material,  has  been  carried 
out.  Some  characteristic  behavior  of  damages  in  the  particulate  composite  materials 
has  been  revealed. 
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